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Abstract 

We study the coupled surface and grain boundary motion in a bicrystal in the context 
of the "quarter loop" geometry. Two types of physics motions are involved in this model: 
motion by mean curvature and motion by surface diffusion. The goal is finding a formulation 
that can describe the coupled motion and has good numerical behavior when discretized. Two 
formulations are proposed in this paper. One of them is given by a mixed order parabolic system 
and the other is given by Partial Differential Algebraic Equations. The parabolic formulation 
constitutes several parabolic equations which model the two normal direction motions separately. 
The performance of this formulation is good for a short time simulation, ft performs even better 
by adding an extra term to adjust the tangential velocity of grid points. The PDAE formulation 
preserves the scaled arc length property and performs much better with no need to add an 
adjusting term. Both formulations are proven to be well-posed in a simpler setting and are 
solved by finite difference methods. 



Key words: Grain Boundary; Mean Curvature Motion; Surface Diffusion; Well-posedness; Finite 
Difference 

1 Introduction 

Coupled surface and grain boundary motion is an important phenomenon controlling the grain 
growth in materials processing and synthesis. A commonly used model to study this coupled effect 
is called "quarter loop" geometry introduced by Dunn et al. [6]. 

In the quarter loop geometry, there are two grains between which there is an interface called 
grain boundary as shown in Fig.l. The two grains are of the same material and differ only in their 
relative crystalline orientation. The grain boundary runs parallel to a free surface before it turns 
up and attaches to upper surfaces at a groove root. When heated at a specific temperature, the 
grain boundary migrates to reduce the surface energy and to heal the orientation mismatch. Since 
the driving force is constant the grain boundary moves at a constant velocity after a short time 
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Figure 1: The quarter loop bicrystal geometry. 



decay. It is reasonable to assume that the bicrystal is uniform along the cross-section direction. 
Thus it is reasonable for us to consider only two dimensional(2D) geometry in this paper. 

This geometry contains two types of motion. One of them is mean curvature motion for the 
grain boundary. And the other one is surface diffusion for the upper surfaces. More detail about 
this model is give in [10]. 

Motion by mean curvature is an evolution law in which the normal velocity of an interface is 
proportional to its mean curvature. More precisely, the motion of an interface T satisfies 

V c = Ak (1) 

Here V c denotes the velocity in the normal direction of T, and k stands for the mean curvature of 

r. 

First proposed by Mullins [13] to model the curvature driven diffusion on the surface of a 
crystal, surface diffusion is a different evolution law in which the normal velocity of an interface is 
proportional to the surface Laplacian of mean curvature. The motion of interface T satisfies 

V d = -BA s k (2) 

Here Vd stands for the normal velocity and A s stands for the operator of surface Laplacian which 
is defined as 

A s = V s ■ V s where V s = V — nd n (3) 

In two dimensions, surface diffusion can be reduced to a normal direction motion with a speed 
function depending on the second derivative of the curvature with respect to arc length, i.e., 

V d = -Bn ss (4) 
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upper surfaces 




Figure 2: An example when one of the surfaces is not single- valued. 



Here s is arc length parametrization. Since the problem is proposed in two dimensions, the motion 
by surface diffusion always refers to equation (4) instead of the general case (2) in this paper. 

We shall prove later in the appendix that we could normalize A and B by rescaling the time 
and space. Since this reformulation make the problem neither harder or easier we will take both A 
and B as one. 

Surface diffusion is an intrinsically difficult problem to solve numerically even in two dimensions. 
The main difficulty is that it is stiff due to the fourth order derivatives and such that an explicit 
time stepping strategy requires very small time steps. Moreover, owing to the lack of a maximum 
principle, an embedded curve may not stay embedded, in other words, it may become self-intersected 
during the evolution. 

Travelling wave solutions have been derived for the whole nonlinear problem and for a linearized 
problem in [11] and [10] respectively. These analytic results are used to verify numerical results in 
this paper. 

The formulations in this paper will be proposed in parameterized form. There are two reasons 
why we prefer the parameterized form. Firstly, we try to set up a versatile formulation that is 
extensible to other problems, such as the evolution of a closed curve , other freely positioned triple 
junction problems and even the evolution curve networks, wherever for a single type motion or a 
mixed type motion. Secondly, even for the coupled grain boundary motion, the function y(x) which 
represents the exterior surface may not be single-valued as shown in Fig.2 and this phenomena is 
physically reasonable [8, 9]. Also for such consideration we will treat the exterior surface as two 
curves separated by the triple junction in the following discussion. 

The outline of this paper is as follows. In section 3, parabolic equations are derived for the 
motion by mean curvature and the motion by surface diffusion separately. The boundary conditions 
including the triple junction conditions and domain boundary conditions are discussed in section 4 
referring to the analytic work of Novick-Cohen et al. [10, 11] and Wong et al. [12]. In section 5, the 
well-posedness for a linear parabolic system that is closely related to the full nonlinear problem is 
analyzed and followed by a discussion about the artificial tangential condition. From section 6 to 8 
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we discuss the numerical details including discretization, time stepping and some other numerical 
issues. We start the discussion for a PDAE system in section 9. An interesting computational 
example of surface diffusion is given in section 10. The linear well-posedness of the PDAE system 
is analyzed in section 11. 



2 Cartesian Formulation 



In this section, we consider the problem in the cartesian coordinate system which is give as be- 
low(see [10]). 
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Here 2/ = 2/( a; ) stands for the height of the two exterior surfaces, u = u(x, t) is the height of the 
grain boundary and s(t) denotes the location of the junction where the three surfaces meet. 

Since the junction is moving it is not straightforward to solve this system numerically. We fix 
the junction by making the following transform, 



x = x — s(t) 

We let y(x,t) = y(x,t) and u(x,t) = u(x,t). Therefor, 

y x {x,t) = y x (x,t) 

y t {x,t) = y t (x,t) - y x {x,t)s t 

u x (x,t) = u x (x,t) 

u t (x,t) = ut(x,t) - u x (x,t)s t 



(7) 
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Result at time 0.1 
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Figure 3: Numerical result for system (8)-(9) with m = 0.5. Dotted line: numerical result; Solid 
line: travelling wave solution. 



Then system (5) becomes 
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And the boundary conditions (6) become 

y(0 + ,t) = y(0-,t)=u(0 + ,t), t>0 
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This system could easily be discretized using standard finite difference schemes on a fixed 
staggered grid. A numerical result is shown in Fig. 3. 
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The disadvantage of using cartesian formulation is that it is not applicable to non-single valued 
case as shown in Fig. 2. What's more, since the grain boundary is nearly singular at the junction, 
it requires very small grid size for accuracy. For wider application we consider two parametric 
formulations in the rest of this paper. 

3 A Parabolic Formulation 

In this section we derive a parabolic system to model the coupled motion. Here and throughout 
this paper we use X = (u(-),v(-)) to represent a parameterized curve with u(-) and v(-) being the 
coordinates. 

Several more notations should be introduced as well. The arc length parameter is denoted 
by X(s) = X(u(s),v(s)) and any other parameter is denoted by X(a) = X(u(a),v(a)). t and ft 
stand for unit tangential direction and unit normal direction respectively, k stands for curvature. 
Although all final equations are parameterized by a, the arc length parametrization is useful for 
the intermediate deviations. 

3.1 Motion by Mean Curvature 

We first derive a parabolic equation to describe the motion by mean curvature. Similar discussion 
has been addressed in [3] and [7]. We give a brief description for reader's convenience. 
With the notations introduced above one has 

X s = t 



X ss = nn 

Here the subscript s stands for the derivative of X with respect to arc length s. Direct computation 
shows that 

x a = x s ^P- = X s \X a \ (10) 
da 

where S(a) is defined by 

S(v) = T y/^ + rtdv (11) 

J CO 

which stands for the length of the curve from point X(gq) to X(a). \X a \ is L2 norm of X a defined 
by 

\X a \ = ^ul + vl 
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Differentiate equation (10) with respect to a to botain 

Xacr = X ss | X a \ 2 + A S |X CT | s | X a | (12) 

By previous derivations one can compute the normal component of vector and obtains 



X aa _ X ss \X a \ s X ss 

n = X S g h -r—rXg ■ 



\X a \ 2 K \X C \ K 

= K (13) 

Thus, if we set up a formulation: 

^ = DTI* ( 14 ) 

it is obvious that the motion described by (14) has normal velocity k. This gives us an option to 
describe the motion by mean curvature. Equation (14) is fully parabolic which means it is parabolic 
in both the normal component and tangential component. 

There are some other equations that can also describe curvature motion, for example, 

X t = Kft (15) 

But this system is not fully parabolic. A linearization shows that this system is parabolic in the 
normal component and hyperbolic in the tangential component. A discretization of system (15) 
will not have the good numerical properties as those of a fully parabolic system due to the lack of 
regularity in the parametrization as shown in [3]. 

3.2 Motion by Surface Diffusion 

Considering the good properties of a parabolic formulation, we hope to find a parabolic formulation 
for motion by surface diffusion. By analogy with the approach to the mean curvature motion 
described above, we try the following form: 

X 

Xt = — | , 4 + L(X aaa , X aa , X a ) (16) 
|^-<t I 

where L(X aaa ,X aa ,X a ) includes some lower order terms and will be determined such that 

X 

x t-n = ( — ; + L(X aaa ,X aa ,X a )) • ft = -K S g (17) 

We focus our study on finding out L(X aaa , X aa ,X a ) in the rest of this section. 
Note first the following equation (see appendix for proof), 

(Xg S gg + K 2 X SS ) ■ n = K ss (18) 
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Compare equation (17) and (18) to get a choice for L, 

L(X aaa , X aa , X a ) = , , 4 — X ssss — k X ss (19) 
1^*7 1 

One will find later that the fourth order terms appeared in (19) could be cancelled with each other 
and such that L involves only third or lower order derivatives. 

Start by equation (10) and differentiate several times with respect to a to get following relations, 

Xaa = X SS S%. + X s S a(7 (20) 

(21) 
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I (J | D (T D cr D cr D <7 

Substitute equation (23) into (19), rewrite arc length parametrization s into cr using (20)-(21). 
Since vector X s is perpendicular to n and has no contribution to normal direction we can ignore 
all X s terms and obtain 

TlV V V \ _ a^ aa Xcracr , r^aa X a a .S aa a X aa 2 X aa 

L\A acc , A CTCT , A CT J — D-^-|^ p - |Jf |2 + ~53~|X | 2 ~~ |X | 2 

Substitute equation (24) back into (16) and collect to get the scheme as 

Y Y <? 2 <? Y 

We would like to point out that the choice of L is not unique. A similar expression has been 
given by Garcke et al. in [7]. 



3.3 The Parabolic System 

We now give the fully parabolic system, 
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Figure 4: Sketch of the grain boundary groove. 



where X 1 stands for the grain boundary and X 2 ,X 5 stand for the left branch and right branch of 
the upper surface respectively. All curves are represented by X(o) with a G [0,oo). 

This system will be solved numerically with the boundary conditions discussed in the next 
section. 

4 Boundary Conditions 

The grain boundary and the two upper surfaces meet together at one end which is referred as 
triple junction. The other end of the three curves tends to infinity in the quarter loop geometry. 
For numerical reasons, we compute this problem in a bounded domain. This domain is chosen 
large enough such that it can simulate the motion at least for a short time. This restriction is 
reasonable since the curves are asymptotically flat for the parts far away from the triple junction. 
All computations presented in this paper are constrained in a finite domain [—6, 12] and the curves 
are parameterized with a G [0, 1]. 

At a = the three curves meet at a triple junction and at a = 1 the three curves meet the 
artificial domain boundary separately. 

4.1 Triple Junction Conditions at a = 

We first discuss the boundary conditions at the triple junction. First of all, three curves should 
have common coordinates at a = 0, i.e., 

X\Q, t) = X 2 (0, t) = A 3 (0, t) (27) 
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By Young's law we have two more conditions which are referred as angle conditions: 



Xl Xl „ . m 



cos 6*i2 = cos( h arcsin — ) (28) 



\Xl\ \Xl\ v 2 2 



xl xl 



ir m 

| X1 | | X 3| cos 6»i3 = cos(- + arcsmy) (29) 

where Oij denotes the angle between curve i, j and m = ^ grain 1 1 exterior is a constant measuring the 
relative surface energy between the grain boundary and exterior surface. 
The continuity of the surface chemical potentials implies that 

* 2 = ~ k3 {k= wj*~ ra } (30) 

Here the superscripts are indices of curves. 
And the balance of mass flux implies that 

«2 = 4 («. = x ™'*° - 3 |XctIct [^; Xct±) ) (3i) 

where the expression for k s is obtained by taking the derivative of the expression of k directly. 

We must be careful about condition (30). Basically, we need the two upper surfaces have the 
same convexity. Since a has opposite directions for the two curves the odd time derivatives will 
have opposite signs when computed by parametric form. Thus we should put a minus sign for (30) 
and keep the same for (31). 



4.2 Boundary Conditions at a = 1 

At the other ends of the curves we put several artificial conditions such that they do not move 
during evolution and keep being flat. This is reasonable since they start being flat and they will 
not be influenced by the motion of the triple junction in a short time. The following conditions are 
imposed at a = 1, 

X t i (l,t)=0 fori = 1,2,3 
Xi a (l,t)=0 fori = 2,3 



4.3 Artificial Tangential Conditions 

We should point out that the whole system contains two second order equations and four fourth 
order equations and it should have ten conditions at the junction point for well-posedness. Recall 
that there are only eight junction conditions as have been addressed above. Thus, we need two 
more conditions. There are several options to impose the extra conditions. And we will prove later 
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that different conditions could only change the parametrization of the curves and will not change 
the profiles of the curves. Since these conditions do change the tangential velocities of the grid 
nodes we refer them as artificial tangential conditions. As one of the options, the following two 
conditions are applied into the system: 



5 Well-posedness for the Parabolic System 

In this section, we analyze the well-posedness of the system proposed above. We linearize around 
fixed straight line solutions and get a system that has the same highest order parabolic behavior 
as the original problem. The well-posedness we do gives the conditions that match those that in 
more complicated nonlinear analysis gives, where such analysis exists. And therefore, we believe 
the results of the analysis should apply to the full nonlinear problem. 

5.1 Linearization of the System 

To linearize the system we consider a perturbation expansion around the tangential direction at 
the triple junction for each curve, i.e., 



where di = (dn,di 2 ) is a constant vector standing for the unit tangential direction. Substitute 
above equations into (26), linearize and keep the leading order terms to get a linear system: 
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For convenience, we omit the bar above X in following discussion. 

The linearization of the triple junction conditions is straightforward. 



• Common point at a = 0: 



X 1 = X 2 = X 3 



• Angle conditions: 



di • X 2 + d 2 ■ xl - (d! • d 2 )(d 1 ■X 1 a + d 2 -X 2 ) = 
di • X 3 + d 3 ■ Xl - (di • daXdj • Xl + d 3 ■ Xl) = 



11 



Continuity of surface chemical potentials: 



• Balance of mass flux: 



Xi 



• Artificial tangential conditions: 
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The linear system (33) can be solved using Laplace transforms to get 
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(34) 
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and here s temporally stands for the transformed time variable of Laplace transform. 

For simplicity, we first suppose the angles between any two curves are |7r. Substitul 
(34) into boundary conditions one obtains a 10 x 10 coefficient matrix M(transposed) 
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Linear well-posedness requires that the determinant of matrix M is nonsingular for any s sat- 
isfying Re(s) > 0. Since the well-posedness depends only on their relative positions, we suppose 
further that 




With these assumptions, one obtains the determinant of M: 

\M\ = 6V6s 11/4 + 24v / 3s 3 

Similarly the determinant of M for arbitrary angles is 

\M\ = 32(sin6»i 3 sin 2 6 12 + sin6»i 2 sin 2 6»i 3 )s 3 - 16>/2(siii 9 12 sm9 13 sm(9 12 + 6>i 3 ))s 11/4 

where 6*12, #13 are the angles between the curves as shown in Fig.4. M is nonsingular for any s 
with Re(s) > if < #12, #13 < 7r. The constraint on 9 is not an issue since it has included all the 
cases of interest. 

5.2 Analysis of Artificial Tangential Conditions 

As have been mentioned before, there are several options for the artificial tangential conditions. We 
are interested to know if different choices will lead to the same solution which is shown to be true. 
To prove this point, it suffices to prove that the position of the junction and the three tangential 
directions do not depend on the artificial tangential conditions. The idea to prove this point is to 
show that they all lead to the same solution for X\. If this is true, the position of the junction 
point and the tangential direction of X\ are uniquely determined. Since the angle conditions are 
guaranteed, the tangential directions of the other two curves could also be uniquely determined. 
To sum up, the key point is proving coefficients of X 1 , i.e., An,Ai 2 do not depend on the extra 
conditions. 

The coefficients Ay,Bij in solution (34) can be solved by 

M-C = P (36) 

where M is the coefficient matrix (35) for boundary conditions, C = [An, A\ 2 , ■ ■ ■ ,^32,-632] is 
the coefficient vector to be solved and P = \pi,P2, ■ ■ ■ ,P9,Pio] is a constant vector depending on 
the initial data. Note that only pq,pio and the last two lines of M depend on artificial tangential 
conditions. 

According to the discussion above we need to prove An,A\ 2 do not depend on the artificial 
tangential conditions. More precisely, we need to prove An,Ai 2 do not depend on the last two 
lines of matrix M and P9,pio- 
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For convenience, we rewrite M into a block form 



Mi (8 x 2) M 2 (8 x 8) 
M = | ! (37) 

M 3 (2 x 2) M 4 (2 x 8) 

We do the Gauss elimination for block M 2 and it shows that the rank of submatrix M 2 is 6 for any 
angle conditions, this means we can make the last two lines of M 2 be zeros by row deduction and 
meanwhile making the last two lines of Mi into a full rank (2 x 2) matrix. 

We again use M to denote the new matrix after row deduction. Next we compute M _1 in a 
block form satisfying 



M x M" 1 



Mi (8 x 2) M 2 (8 x 8) \ / Mi (2 x 8) M 2 (2 x 2) 

M 3 (2 x 2) M 4 (2 x 8) / \ M 3 (8 x 8) M 4 (8 x 2) 
1(8 x 8) 

1(2 x 2) 



(38) 



Expand directly to get 



Mi x Mi + M 2 x M 3 = 1(8 x 8) (39) 
Mi x M 2 + M 2 x M 4 = 0(8 x 2) (40) 

Note that (39)- (40) do not involve M 3 ,M 4 which means they do not depend on the artificial 
tangential conditions. If Mi, M 2 can be determined by equation (39)-(40) then we can say Mi, M 2 
do not depend on the artificial conditions. The fact 



An 

Ml 



\ = (M 1 M 2 )xP 



implies that An, Ayi do not depend on the artificial conditions if we can further prove M 2 = . 

Actually, Mi can surely be solved from equation (39). This is because the last two lines of M 2 
are zeros and we have exactly sixteen equations involving only the sixteen unknowns of M\. For 
the same reason we can solve for M 2 by equation (40). Actually, since the last two lines of Mi is a 
full rank (2 x 2) matrix M 2 must be zero. This completes the proof that the coefficients An, A\i in 
equation (33) do not depend on the artificial conditions. And consequently, the shapes of the three 
curves do not depend on the artificial tangential conditions. Novick-Cohen et al. [9] also pointed 
out that the artificial conditions do not influence the solutions, although the problem there is a 
little bit different. In [9] the authors look at a three phase problem in which all three interfaces 
evolve by minus the surface Laplacian of mean curvature and meet at a triple junction. 
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6 Numerical Discretization 



Back to the full nonlinear problem, we present in detail the discretization procedure of the parabolic 
scheme (26) and junction conditions (27)-(31). The basic approach is to use a staggered grid in a and 
we shall denote the approximations by capital letters with subscripts, i.e., Xj(t) ~ X((j— l/2)h, t) = 
(u((j — l/2)h,t),v((j — l/2)h,t)) where h is grid spacing and N = 1/h is the number of interior 
grid points for a £ [0, 1]. 

In order to write the discretized equations we introduce some additional notations. Let 
denote the second order centered approximation of the kth derivative, i.e., 

D 1 X J = (X j+1 -Xj-{)/2h 
D 2 Xj = {X j+1 + X j _ 1 -2X j )/h 2 

and let D + and T denote forward differencing and forward averaging, respectively, 

D+Xj = {X j+1 -Xj)/h 
TX 3 = (X j+1 + Xj)/2 

We discretize each motion separately. 

6.1 Grain Boundary Motion(Motion by Mean Curvature) 

The grain boundary motion is approximated at all grid points by standard differences, 

D 2 X) 

X\ = - 1- i = l,j = 1,2,- ■■ ,N (41) 

j |D 1 X*| 2 

where Xj stands for time derivative. Formally, these discrete equations require values of Xq and 
Xjv + i outside the computation domain. We shall use the boundary condition to extrapolate the 
interior values of X\ and X^ to the unknown exterior values of Xq and Xjy+i. We shall give the 
details of the extrapolation procedure later. 

6.2 Surface Diffusion 

The higher order derivatives appeared in surface diffusion are approximated by 

(X_), . D 3Xj = ^X 3+l -D 2 X 3 ^ (42) 



(X^) 3 c D 4 X, = ° 2X ^ + D f 2 ]+1 ~ 2D2Xj ( 
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triple junction 




•grain boundary 



Figure 5: Sketch of the ghost points at the triple junction. 



There are some other terms such as S a , S aa , S aacT to be approximated. Start from (11) and differ- 
entiate several times with respect to a to get 



S a = yju 2 a + vl = \X a \ 

^a^aa ^a^aa X a ' X a a 



[X c ■ Xaa) 2 X aa ■ X aa + X a • X aatT 



\X I 3 \X I 

I & I I c| 

Every term in scheme (26) is now ready to be approximated by standard differences. 
6.3 Junction Conditions at a = 

The discretization at the junction point is much more complicated. Since there are fourth order 
derivatives for the surface diffusion we shall need two ghost points for each surface curve and 
one ghost point for grain boundary. These ghost points are denoted by Xq , X^_ ± , Xq , X 3 _ 1 , Xq 
respectively. The junction conditions (27)-(31) are approximated as follows, see Fig. 5. 
Condition (27), 

FXl = TXl = TXl = C (44) 

where C denotes the junction point. 

The angle conditions (28)-(29) are approximated by 

D+Xq D+X% . . 

" -= cos 6i2 (45) 



\D + Xl\ \D + X%\ 
D+Xl D + X$ = cog0i3 



\D+X*\ \D + X$\ 
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Discretize condition (30) for each surface curve to get 

(47) 



D 2 X 2 C ■ (DiX^y- _ £> 2 *% • Pi*c) 



Since staggered grid are used, center C is a midpoint not a grid points. But we still can use previous 
notations Dk with the following extensions 



k s can be expressed by 



Xb-i = {Xl 1 +X i Q )/2 = TXU 
Xh +l = {X\+Xi)/2 = FX\ 

Ks = Xa ™'*° - 3 5 ^ ( ^ 5 ^ } (48) 



Thus condition (31) is approximated by 

D 3 X* ■ (DjX*^ JUD2X 2 C -{DiXl)^) = D 3 Xl,.(D 1 X^ 5g g (£> 2 Xg, . (D.X^) 
\ D i x c\ 4 \ D i x c\ 5 l^i^cl 4 " l^i^cl 5 

Finally, the artificial tangential conditions is calculated by 

D 2 X}, • DiXh = /oH = 2 3 (5Q) 

We now finish discretizing the junction conditions. 

6.4 Domain Boundary Conditions at a — 1 

The discretization at <r = 1 is straightforward. 

•^ArLn-d* = jrX ^lt=(n-i).dt for * = X > 2 ' 3 
D+X^ = fori = 2,3 

7 Time Stepping 
7.1 Explicit Scheme 

As an explicit scheme, forward Euler method is used for the time stepping process. 

x n+l =X n + AtF ( X n ) (51) 

Here F(X n ) denotes the right hand side in formulation (25) evaluated at time level n. Time steps 
At are chosen so that the full discrete scheme is stable. Here we choose At = le — 12. This scheme 
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is easy to implement. Given the results at time n we update the values of the interior grid points 
by forward Euler method to time level n + 1 for the three curves respectively. Next solve for the 
ghost points, junction point and the boundary points by the boundary conditions. Then go on 
to the next time level. The time step is excessively small due to the stiffness of the fourth order 
parabolicity as noted previously. 

7.2 Implicit Scheme 

In order to avoid the excessively small time steps for explicit scheme we consider implicit techniques 
in this section. For simplicity we use backward Euler method. Given the values at time n we update 
the values at time level n + 1 by solving the nonlinear system 



Since the three curves are strongly coupled by the junction, we solve all unknown points simul- 
taneously including the ghost points and the extrapolated boundary points. This leads to a large 
nonlinear system which is solved by Newton's method. There is no doubt that this scheme should 
be stable for any time steps. But it can not survive a long time computation due to the nonuniform 
tangential velocity which leads to a nonuniform distribution of the grid points. This phenomenon 
can not be fixed even if we refine the grid. 

One way to overcome this difficulty is regridding the grid points once they become too far or 
too close. But the bad distribution could happen only near the junction and the closer to the 
junction the sparser (or denser) the grid points are. Hence it is hard to regrid no matter globally 
or locally. Another way is adjusting the tangential velocity of the grid points such that they could 
adjust themselves being uniform. And this is the motivation for the next section. 

Numerical results for scheme (26) with time step At = le — 4 are shown in Fig. 6. All numer- 
ical experiments in this paper start from the same position as showed in Fig. 6. All results arc 
compared with a travelling wave solution solved by Amy et al [11]. 



8 Adjustment of Tangential Velocity 

We have mentioned that long time computations are problematical even for implicit schemes. This 
is because of the bad distribution of grid points. To get a more uniform distribution of gird points 
along the curve we consider adding an artificial term to adjust the tangential velocity of the grid 
points for the motion by surface diffusion. 

We consider the following modified scheme for the fourth order problem 



X n+1 - AtF(X n+1 ) -X n = 



(52) 



X t = F{X) + a{ 




(53) 
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initial status Result at time 0.2 




-6 -4 -2 2 4 6 8 10 12 -3 -2 -1 1 2 3 4 



Figure 6: Plot of results for scheme (26) with backward Euler method for a short time with 
m = 0.5. Left: initial status with grid points; Right: result zoomed in near triple junction. Dotted 
line: numerical solution; Solid line: travelling wave solution; Time step size: At = le — 2. 

where a is a constant to be determined. The newly added term in (53) will not influence the 
normal velocity but it does change the tangential velocity. We do not know exactly how to choose 
the optimal a but a = —100 seems to work well for our problem. The result is shown in Fig. 7. It 
is obvious that the grid points are much more uniform than that in Fig. 6. The time step size for 
Fig. 7 is At = 0.01. A numerical convergence study is shown in Table 1. 



dt 


As 


L2 Norm 


Rate 


Lqo Norm 


Rate 


dt = O.OlAs 2 


0.2 


3.1494e-04 




2.0241e-03 




0.1 


7.9775e-05 


1.9811 


5.4797e-04 


1.8852 


0.05 


2.1445e-05 


1.8953 


1.4530e-05 


1.9151 



Table 1: Estimated errors and convergence rates for parabolic formulation with m = 0.5. Errors 
are evaluated at t = 0.02 

Although the newly added tangential term improves the numerical behavior it can not com- 
pletely overcome the difficulty The artificial tangential conditions discussed in section 5 make the 
problem even more complicated. All these motivate us to seek a more efficient scheme. 
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Result at time 0.2 Result at time 0.2 




Figure 7: Plot of the results for scheme (53) with m = 0.5, a = —100, At = 0.01. Left: result at 
t=0.2; Right: result zoomed in near triple junction at t=0.2. Dotted line: numerical result; Solid 
line: travelling wave solution. 

Remark There is another way to adjust the tangential velocity, 

I f = f(I) + «i.t)t (54) 

\-X-a\ 

For this case, we should choose a positive, for example a = 100. 

9 A PDAE Formulation 

As we have pointed out in section 8, the fully parabolic scheme (26) does not always have good 
numerical behavior. And the presence of the artificial tangential condition makes the discretization 
of the original problem more complicated. In this section we propose another formulation that 
can overcome these disadvantages and also avoids possible loss of tangential monotonicity in the 
parametrization due to the fourth order PDE. 

Let us again start from the motion by mean curvature. First of all the basic evolution law 
should be satisfied, i.e., 

X t ■ n - k = (55) 

Because there are two free variables in this equation we need one more equation for solvability. 
Since the requirement for the normal direction motion has been fulfilled by equation (55) we use 
the second equation to impose a constraint on the distribution of grid points. It is natural to let 
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all grid points have equal spaces. To avoid introducing an extra variable we let the change rate 
between any two adjacent spaces is zero, i.e., 



(56) 
Note that 

\X a \ a = (\/ X a ■ X a ) a = — — 

\^a\ 

the following equations are actually used to describe the motion and keep grids equi-spaced, 

Xt • n — k = 
Xa ■ X cc = 

These are called partial differential algebraic equations (PDAEs). 

In a similar way we derive the PDAEs for the motion by surface diffusion, 

X t ■ n + k ss = 
X c ■ X aa = 

Then the full PDAE system for the coupled motion is 

X]-n-K = Q 
Xl ■ n + k ss = 

A t 3 • n + k ss = (57) 
Xi-Xi a = i = 1,2,3 

The boundary conditions are the same as the parabolic case except that we do not need artificial 
tangential conditions any more. 

This is an implicit index- 1 DAE system. Usually an index- 1 DAE can be discretized directly 
without any numerical difficulties [1], and that is our experience in this case. 

Although the boundary conditions are the same as those of the parabolic system, the discretiza- 
tion is a little bit different. Instead of using five ghost points we now introduce only three ghost 
points plus two extra variables which stand for the curvature at the two ghost points corresponding 
to the two surface branches. The two variables are denoted by Kq,Kq and the last two junction 
conditions are approximated by 

' 1 ' 1 

2 ~ 2 
(kq — ft|) (kq — nf) 

\D.XH 
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where k\ stands for the curvature of the first interior point of curve i and we use the average of 
/cq, k\ to approximate the curvature at the center point, i.e., junction point. Again the sign should 
be carefully handled. 

Implementing this scheme one obtains a better result shown in Fig. 8. The result is much more 
accurate and the grid points are more uniform as well. An error comparison is shown in Table 2. A 
numerical convergence study of the PDAE formulation is shown in Table 3. The convergence rates 
shown in Table 3 are close to 2 as expected. 



Result at time 0.5 Result at time 0.5 




-1 -0.5 0.5 1 1.5 -1 -0.5 0.5 1 1.5 



Figure 8: Results comparison between the two schemes. Left: result for (26); Right: result for 
(57). Both pictures are zoomed in near triple junction. Dotted line: numerical solution; Solid line: 
travelling wave solution; Time step size: At = 0.01. 





Parabolic Formulation 


PDAE Formulation 


As 


0.05 


0.05 


At 


0.01 


0.01 




0.0041 


0.0027 



Table 2: Performance of the two formulations with L ro norm and At = 0.01. 

Without difficulty we can apply this scheme to the case when the surface curve is not a single- 
valued function as shown in Fig. 9 
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dt 


As 


L2 Norm 


Rate 


Lqq Norm 


Rate 


dt = O.OlAs 2 


0.2 


2.7837e-04 




1.8996e-03 




0.1 


7.2717e-05 


1.9366 


5.4444e-04 


1.8029 


0.05 


1.8732e-05 


1.9568 


1.4732e-04 


1.8858 



Table 3: Estimated errors and convergence rates for PDAE formulation with m = 0.5. Errors are 
evaluated at t = 0.02. 



Result at time 0.2 




-0.5 0.5 1 1.5 2 



Figure 9: Plot of the results for scheme (57) with m = 1.96 which has non-single valued upper 
surface. Dotted line: numerical result; Solid line: travelling wave solution. 

10 An Example of Surface Diffusion Problem 

We temporally move our focus to a normal direction motion that involves only motion by surface 
diffusion. The motion starts with a closed star shaped curve and evolves with a speed equal to the 
second derivative of curvature with respect to arc length. According to the properties of surface 
diffusion the curve will evolve into a circle and preserve the area enclosed by itself. This problem is 
computed using the PDAE formulation for the surface diffusion and the result is shown in Fig. 10. 
The method conserves the area quite well and the change is about 0.032%. Similar examples have 
been investigated using level set methods in [5, 14]. Level set methods have unbeatable superiority 
for interface motion problem especially when there is topology change. But for this simple problem 
(with no topology change ) our method is more efficient and accurate. Note that this problem can 
also be computed by the parabolic formulation. 
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Time 



Time 7.5e-005 





Time 0.0001745 



Time 0.0003575 
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Figure 10: A computational example that involve only the motion by surface diffusion. At 
5 x 1(T 7 . The area changes by 0.032% 



11 Well-posedness for the PDAE System 

Similar to the parabolic system we do a well-posedness analysis for the PDAE system in this section. 
Considering the same linear problem as that in section 5.1 one obtains 



x\ 


■di 


— xL ■ di 


di- 


xL 


= 


x? 


•4 


= ~ X Laa ■ 4 


d 2 - 


xL 


= 


x? 


■di 


- Y 3 rl 1 

— ~ JV aaaa °3 


d 3 - 


xL 


= 



where di and df- stand for unit tangential direction and unit normal direction of the i th curve 
respectively. 
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(59) 



Linearization of the boundary conditions are exactly the same as before. They differ only for 
the discretization procedure. 

If dn,di2 7^ the linearized system (58) has solution in the form 

ui = A n e-^ + B u 
ui = -hAue-^ + ±B U 
w 2 = A 2 ie AlCT + B 21 e X2(T + C 2 i 
^2 = -h(A 22 e x ^ + B 22 e A2C7 ) + ±C 21 
u 3 = A 31 e x ^ + B 31 e x * a + C 31 
{ v 3 = -k 3 {A 32 e x ^ + B 32 e x ^) + ±-C 3l 

where h = ^r- is a constant and 

\i = (-— + —^s X 2 = ( - 

Without changing the well-posedness of the problem we specify one of the tangential directions, 
say d\ = (0, — 1) T . Further we assume 

7T 



012,013 G (M) and ^12,013 / j 



Since dn = now the solution is changed to 



vi = Bu 

u 2 = A 21 e XltT + B 21 e x * c + C 21 

v 2 = -k 2 (A 22 e x ^ + B 22 e x ^) + ±C 21 

u 3 = A 31 e x ^ + B 3l e x *° + C 31 

v 3 = -k 3 (A 32 e x ^ + B 32 e x ^) + ±C 31 



(60) 



Apply these solutions to boundary conditions to get an 8x8 matrix M and compute the determinant 
of M directly to get 



\M\ = 



4\/2s 7 / 4 sin(0i 2 + 0i 3 ) - 8s 2 (sin6> 12 + sin 6 13 ) 



COS 2 012 COS 2 013 

For the special case when one of the angles 0i 2 ,0i3 is ^, for example, 0i2 = f , 

( m = A n e-^ a 
vi = Bu 
u 2 = C 2 i 

v 2 = A 22 e Xia + B 22 e X2a 
u 3 = A 3l e x ^ + B 31 e x * a + C31 
{ v 3 = -k 3 {A 32 e x ^ + B 32 e x ^) + ^C 3 i 



(61) 
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The determinant of the coefficient matrix is 

lwl -8s 2 sin6li3 + 4^/2s 7 / 4 cos6li3-8s 2 
M = i2 w ^ 

COS^ 6*13 

And if both #12, #13 are |, 



< 





= Aue-^*" 


Ul 


= B n 


""2 


= C21 


«2 


= A 22 e Al<T + B 22 e X2a 


^3 


= C31 


V3 


= A 32 e Al °- + J B 32 e A2<T 



(62) 



The corresponding determinant of the matrix M is 

|M| = -16s 2 

Similar as before all cases discussed above are well-posed if < #i 2 , #13 < n and #12 + #13 > 7T- Note 
that the well-posedness property here coincides with that for the parabolic system we got before. 

12 Conclusion 

We proposed two formulations to describe the coupled surface and grain boundary motion. Both of 
them are well-posed and easy to be implemented by finite difference method. Numerical results are 
shown to be accurate. The PDAE formulation behaves better than the parabolic form does. And 
since all grid points are equispaced for the PDAE formulation it is convenient to regrid globally 
when necessary. This often happens when the curves expand or shrink quickly. 

It is obvious that these schemes can also be used to simulate the motion of a curve that involves 
only mean curvature motion or surface diffusion as shown in section 10. And they are extensible 
to any normal direction motion. Wherever applicable these methods are more efficient comparing 
to the level set methods. But they can not manage topology changes during the evolution. 
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Appendix A: Reformulation of the Motion. 

The original curvature motion and surface diffusion are given as 

V c = Ak 
V d = -Bk ss 

Without loss of generality will shall prove in particular in a parameterized form that we may 
normalize both A and B by rescaling the time and space. 
Suppose the original motions are modelled by 

Xt ■ n = Ak, 

Y t -n = -Bk ss (63) 



We first rescale the time by 



Then equation (63) becomes 



t = Tt 



Xz • n = —k 

L rp 



Y r n = ~k m (64) 



The spatial variable X, Y are rescaled as 

X = RX, Y = RY 

And we can easily verify that 

s = Rs 

where s is the arc length in the new system. More derivation shows that 



Xg — X s , Y§ — Y s , R — ^> K ' — R3^ ss 



Now equation (64) becomes 



Y - AR2 ~ 

Y r n = — t^Kss (65) 
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Here we use the same notation n since the normal direction does not change. If we choose 



then equation (65) becomes 



V B B 



Xi • n = k 



Y~ t -n = -hss (66) 
This complete the normalization of coefficients A and B. 

Appendix B: Proof of Equation (18) 

One has the following fact 

X ss ■ X ss = k 2 (67) 
Differentiating (67) with respect to s one obtains 

X.sss ' X ss — HLK S 

Xgsss ' X ss -\- X sss • Xggg — fcfcgg ~t~ Kg 

Then an expression for be derived, 

Xgggg • X S g + Xggg • Xggg — k /.„> 
K S g = (68) 

K 

On the other hand, one has 

X ss = nn 

Take derivative again and calculate the inner product of X SS g 

■^sss — K>gTl ~\~ K/Tlg — fcgTl hv t 

Xggg -Xggg = k 2 s + k a (69) 
Substitute equation (69) into (68) to get 

K ss = Xgggg ■ n + k 3 (70) 
Using equation (70) together with the fact that 

k 2 X ss ■ n = K 2 (X SS ■ ft) = k 3 (71) 

one obtains 

Kss = {Xgggg + k X S g) • n 

which completes the proof. 
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